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For rime ice - where the ice buildup has only rough and jagged surfaces but no protruding horns - this study 
shows two-dimensional CFD analysis based on the one-equation Spalart-Almaras (S-A) turbulence model to 
predict accurately the lift, drag, and pressure coefficients up to near the stall angle For glaze ice - where 
the ice buildup has two or more protruding horns near the airfoil’s leading edge - CFD predictions were 
much less satisfactory because of the large separated region produced by the horns even at zero angle of 
attack. This CFD study, based on the WIND and the Fluent codes, assesses the following turbulence 
models by comparing predictions with available experimental data: S-A, standard k-e, shear-stress 

transport, v 2 -^ and differential Reynolds stress. 


1. INTRODUCTION 

When an aircraft flies into environments, where 
meteorological conditions can cause ice to form on its 
wings, the aircraft’s ability to maintain flight will 
diminish quickly with time unless there are ways to 
eliminate the ice formed. This is because with ice 
buildup, which occurs mostly on the wing/airfoil’s 
leading edge, not only is the lift reduced, but also stall 
will occur at much lower angles of attack. 1 If the 
maximum lift force that can be generated by the aircraft is 
less than its weight, then the aircraft will fall from the 
sky. Even if the lift is sufficient to sustain flight, uneven 
ice buildup on the wings can produce flight control 
problems. Thus, it is critically important to understand 
the different ice shapes that can form on the wings and 
how they affect aerodynamics. 

The effects of accrued ice on aircraft aerodynamics 
can be studied by flight tests, wind-tunnel measurements, 
and computational fluid dynamics (CFD) simulations. 
Flight tests are the most realistic but are expensive. Tests 
in wind tunnels offer a controlled environment, but cannot 
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reproduce actual flight conditions since not all 
dimensionless parameters can be matched. CFD is the 
most cost effective method and has the ability to simulate 
actual geometric and flight conditions. However, its 
accuracy depends on the quality of the grid and the ability 
of the turbulence model to reproduce the key flow physics. 

Previous Work 

When studying aerodynamics of wings with ice 
accretion by using CFD, the generation of high-quality 
structured grids is a major challenge. This is true even for 
two-dimensional (2-D) analyses involving 2-D ice shapes 
on 2-D airfoils. Chi, et al 2 presented methods to generate 
high-quality single- and multi-block structured grids for 
complicated 2-D ice shapes. They demonstrated their 
methods on a naturally laminar flow airfoil (NLF0414 1 ) 
with a glaze ice (623 ice shape 1 ), which had two 
protruding horns. To generate high-quality multi-block 
grids, they showed that a “thick wrap-around” grid is 
needed to ensure that grid points clustered next to solid 
surfaces do not propagate into the interior of the flow 
domain (Fig. 1). To minimize and confine the adverse 
effects of the ice’s jagged surface on the quality of the 
grid, they suggested using a transitional layer next to solid 
surfaces (Fig. 2). They also developed a method for 
generating single-block grids for iced airfoil based on an 
idea presented by Tai. 3 

Since multi-block grids can converge slower, Zhu, et 
al. 4 studied the effects of blocking strategies on the 
convergence rate to steady-state solutions. They showed 
that when a single-block grid is partitioned into a multi- 
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block grid with block boundaries perpendicular to the 
flow direction, the number of iterations needed to 
converge is about the same as that for the single block 
grid. When that same single-block grid is partitioned into 
a multi-block grid with block boundaries parallel to the 
flow direction (e.g., wrap-around grids), convergence rate 
can deteriorate if those block boundaries are too close to 
the airfoil/ice surface. If those boundaries are sufficiently 
far away, then the convergence rate can be nearly as good 
as single-block grids. 

Zhu, et al. 5 examined the grid-generation and the 
blocking techniques of Chi, et al. 2 and Zhu, et al. 4 by 
applying them to a much more complicated ice shape, the 
145m ice shape 1 , which has multiple and highly extended 
horns in close proximity to each other. Even for this 
highly challenging ice shape, they found that high-quality 
single- and multi-block grids can be generated but local 
elliptic smoothing is needed (Fig. 3). They also found 
that variable thickness wrap-around grids are needed for 
multi-block grids. On convergence rate, they found that 
one needs to monitor not just the residual for the flow 
variables, but also the residual for the turbulence 
quantities. Though the lift and the drag coefficients 
converge when the flow variables converge, the flow 
patterns were found not to converge until the turbulence 
quantities converged. The residual for the turbulence 
quantities were found to converge slower than the residual 
for the flow variables. 

To date, most CFD studies on 2-D iced airfoils have 
focused on the following: grid generation methods, grid 
resolution, turbulence models, and lift and drag 
coefficients as a function of the angle of attack (see, e.g., 
Refs. 2, 4, 5, and 6 to 1 1). Though much has been learned 
from these studies, accuracy of CFD predictions is still 
unclear. For example, what are the error bounds and 
confidence level in the computed lift and drag coefficients 
as a function of angle of attack? In addition, relatively 
little emphasis has been made on understanding the flow 
field induced by the ice accrued on the airfoil and how 
well they are predicted. One reason is that very little 
experimental data are available that can be used to assess 
the accuracy of such CFD predictions. Recently, 
experimental data for the x-component velocity have been 
made available for a business-jet airfoil (GLC305 1 ) with a 
944 glaze ice shape 1 and is reported in Ref. 12. The 
availability of this experimental data enables a more 
thorough interrogation of the CFD results generated. That 
is, in addition to lift, drag, and pressure coefficients, the 
detailed flow field can also be interrogated with 
confidence to provide better understanding of how ice 
shapes affect aerodynamics. 

Objective and Approach 

With the above backdrop, the main objective of this 
study is to understand how well CFD can predict lift, 
drag, surface pressure, and the velocity field as a function 


of the angle of attack for 2-D iced airfoils. The accuracy 
of CFD predictions will be assessed by comparing 
computed results with experimentally measured lift, drag, 
surface pressure, and velocity field. 

Since it is now possible to generate high-quality grids 
for geometrically complicated 2-D iced airfoils through the 
work described in Refs. 2, 4, and 5, this study on accurate 
CFD predictions focuses on the effects of turbulence 
models, the other major source of error. Of particular 
interest is how well state-of-the-art models can predict the 
aerodynamics induced by glaze and rime ice shapes. 
Glaze and rime ice shapes, formed under different icing 
conditions, constitute the two fundamental ice shapes. 
Though both ice surfaces are rough and jagged, glaze ice 
also has horns but rime ice does not. Thus, these two ice 
shapes produce very different flow fields, which may have 
different requirements on turbulence modeling. For rime 
ice, which produces only very small separated regions at 
all angles of attack except near stall, a simpler turbulence 
model may be adequate. For glaze ice, which can produce 
large separated regions downstream of the horns even at 
zero angle of attack, one- or two-equations models may be 
inadequate. Differential Reynolds stress models that can 
account for each of the Reynolds stresses individually may 
be needed to predict the effects of streamline curvature and 
the time-lagged response of turbulence to changes in the 
mean flow. 

In order to examine as many turbulence models as 
possible, two codes were used in this study. One is the 
widely used, open-source code, referred to as WIND, 
which was used in the work reported in Refs. 2, 4, 5, 7, 8, 
10, and 11. The other is Fluent, a popular commercial 
CFD code with many turbulence models. References on 
WIND and Fluent are provided in Section 3, when these 
codes are described in more detail. 

The organization of the remainder of this paper is as 
follows. Section 2 summarizes the glaze and rime, iced- 
airfoil problems studied. Section 3 outlines the 
formulation and the numerical method of solution used in 
the codes and the turbulence models examined. Section 4 
presents the results generated. 

2. PROBLEM DESCRIPTION 

The airfoil, the ice shapes, and the flow conditions 
selected for study are those for which the flow field is 
sufficiently complicated and for which there are 
experimental data that can be used to assess the accuracy 
of the CFD predictions. The airfoil selected is the 
business-jet airfoil (GLC305 1 ). The glaze ice selected is 
the 944 ice shape with two large protruding horns. The 
rime ice selected is the 212 ice shape, which has 
considerable surface jaggedness but no protruding horns. 
The airfoil and the ice shapes about the airfoil’s leading 
edge are shown in Figs. 4 and 5. See Ref. 1 for details of 
the geometry. 



In this study, the freestream Mach number (M) is 
0.12. Two freestream static pressures (P = 20.5 psi and 
37.0 psi) and two Reynolds numbers based on the free 
stream conditions and the chord length (Re = 3.5 x 10 6 
and 6.0 x 10 6 ) were investigated. The angles of attack 
(AO A) simulated are 0, 4, 6, 7, 8, 9 or the 944 ice and 0, 
4, 6, 7, 8, 9, 10, 11, 12 for the 212 ice. The reason for 
simulating so many angles of attack is to get the details on 
lift and drag coefficients about the stall angle. 

3. FORMULATION AND NUMERICAL METHOD 
OF SOLUTION 

Two different CFD codes were used to generate 
solutions for the iced-airfoil problems described in the 
previous section. One is a widely used, open-source code, 
known as WIND. 13 ' 14 The other is a popular commercial 
code, Fluent. b These codes were selected because they 
are highly versatile and contain a wide range of 
turbulence models. 

For both codes, the flow past the GLC305 airfoil with 
the 944 and 212 ice is modeled by the ensemble-averaged 
conservation equations of mass (continuity), momentum 
(full compressible Navier-Stokes), and energy for a 
thermally and calorically perfect gas. For most 
simulations, the one-equation Spalart-Allmaras (S-A) 
model 16 is used to mimic the effects of turbulence. For 
iced airfoils, Chuang & Addy 10 showed the S-A model to 
out-perform two-equation turbulence models including 
the highly regarded shear-stress transport (SST) 
model. 17 ’ 18 For Fluent, the following turbulence models 
were investigated: S-A, 16 SST, 17,18 standard k-B, 19 

Durbin’s v 2 -f model, 20 and a differential Reynolds stress 
model. 21,22 With Fluent, the near-wall region is always 
modeled by the two-layer model of Chen and Patel. 23 In 
all simulations with both WIND and Fluent, the 
conservation equations and the turbulence models are 
integrated to the wall, where the no-slip and adiabatic 
wall conditions were imposed (i.e., wall functions were 
not used). 

The numerical methods of solution used are as 
follows. For WIND, the convective terms were 
approximated by second-order Roe upwind differencing. 
Since only steady-state solutions are of interest, time 
marching to steady state was accomplished by an implicit 
method based on ADI-type approximate factorization 
with local time stepping. For Fluent, which uses a finite- 
volume method, fluxes at the cell faces are interpolated by 
using second-order upwind differencing. The SIMPLE 
algorithm was used to generate steady- state solutions. 
For this segregated solver, the convergence criteria used 
are to ensure that normalized residual is less than 10' 6 for 
the energy equation and less than 10' 3 for all other 
equations. 

To ensure proper comparison between the codes and 
among the turbulence models, both WIND and Fluent 


used essentially the same grid systems as explained below. 
For WIND, all grid systems generated consist of two 
overlapping single-block grids. One is a fine grid next to 
the airfoil, extending 0.6 chord length from the airfoil in 
all directions (referred to as the inner grid). The other is a 
coarser grid (125 x 21) that overlaps the fine grid by 0.1 
chord length and extends 15 chord lengths away from the 
airfoil in all directions (referred to as outer grid). The 
inner grid is the single-block grid of interest. While 
generating this grid, grid lines were clustered next to the 
airfoil surface so that the first grid point is within a y+ of 
unity. Along the airfoil surface, equal arc-length was 
employed to create a grid as smooth as possible. 

Since Fluent does not accept overlapping grids, the 
grids used by WIND and Fluent are not exactly the same. 
The inner grid used by WIND is also used in Fluent. The 
outer grid used, however, does differ from the one used by 
WIND. In Fluent, the outer grid is an unstructured grid in 
which the square root of the cell area is comparable to the 
grid spacing of the outer grid used in WIND. Since the 
grids closest to the iced airfoils are identical up to 0.6 
chord length for both codes and that is where all of the 
interesting flow features take place, there is a basis to 
compare the two codes. 

In this study, five different inner grids were used. For 
all five inner grids, the same outer grid is used - structured 
for WIND and unstructured for Fluent. The details of the 
inner grids are as follows: For the business-jet airfoil 

without ice (referred to as clean), the inner grid has 913 x 
101 grid points (Fig. 6). The second grid (referred to as 
S25) has the 944 ice shape represented by only 25% of the 
control points (not shown but similar to the one shown in 
Fig. 7 top). This means that there is smoothing of the 
jagged ice geometry, which makes grid generation easier. 
For this grid system, the inner grid has 941 x 101 grid 
points. The third grid (referred to as SV1) uses 100% of 
the control points (Fig. 7 top). For this grid, the 944 ice 
shape was smoothed less so that grid generation is more 
tedious in having to capture more details of the jagged 
geometry. For this grid system, the inner grid also has 941 
x 101 grid points. The fourth grid (referred to as SV2) also 
uses 100% of the control points. SV2 differs from SV1 in 
having more grid points in the region next to solid surfaces 
(989 x 129 grid points). The fifth grid generated is for the 
212 ice, and it has 987 x 131 grid points (Fig. 7 bottom). 

All grid systems described above were generated by 
using transfmite interpolation 24,25 in the manner described 
in Refs. 2, 4, and 5 with varying degrees of local elliptic 
smoothing. Also, all of these grid systems are arrived at 
after an extensive grid-independent study. 

4. RESULTS 

The main objective of this study is to assess how well 
CFD can predict lift, drag, surface pressure, and the 
velocity field as a function of the angle of attack for a 2-D 



airfoil with a glaze or a rime ice shape built upon its 
leading edge. The focus is on the effects of turbulence 
modeling on the predictions. 

Clean Airfoil 

Figures 8 and 9 show the results obtained by using 
the WIND code with the S-A turbulence model for the 
clean GLC305 airfoil. These figures show that CFD can 
predict the lift and surface pressure coefficients quite well 
for angle of attack up to near stall. But, the stall angle is 
slightly under predicted. One reason for not predicting 
the lift coefficient correctly near stall is that the flow 
under those conditions becomes unsteady. In this study, 
only steady-state solutions were sought. 

Airfoil with Glaze Ice 

Figures 10 to 12 show the results obtained by using 
the WIND code with the S-A turbulence model for the 
GLC305 airfoil with the 944 glaze ice. These figures 
show CFD to under predict the lift and drag coefficients at 
all angles of attack. The predicted pressure coefficient is 
shifted towards the leading edge, indicating incorrect 
predictions of the separated region induced by the horns. 

To ensure that grid resolution or ice-shape smoothing 
were not the cause of the inaccuracies, simulations were 
performed with grids based on smoothed (S25) and 
unsmoothed (SV1) ice shapes as well as increased grid 
resolution (SV2). Figure 13 shows these effects to be 
insignificant. Note that SV1 is already a very fine mesh, 
and is the one used in all remaining simulations of the 944 
ice with the WIND and Fluent codes. Thus, using a 
smoothed ice shape (e.g., one represented by only 25% of 
the control points) may be adequate when predicting lift 
and drag. This can be significant, since grid generation is 
easier and less time consuming with a smoother ice 
surface. 

In an attempt to improve predictions, more 
sophisticated turbulence models were evaluated by using 
the Fluent code, which has more turbulence models 
encoded. Figures 14 and 15 show that WIND and Fluent 
gave nearly the same results when the S-A model is used 
on the same inner grid. This gives some confidence 
towards using two different codes to evaluate a variety of 
turbulence models. Figures 14 and 15 show that S-A 
gives the best results for the lift and drag coefficients, 
confirming the findings of Chung & Addy. 10 SST and the 
standard k-e models were found to under predict lift and 
to over predict drag when compared to the other models. 
Differential Reynolds stress and v 2 -f models gave the best 
results for the drag coefficient at AOA = 4° but not at 
AOA - 6°. Also, the lift is severely under predicted at 
AOA = 6°. The less than satisfactoiy results for the v 2 -f 
and the differential Reynolds stress models is that Fluent 
uses the one-equation Chen and Patel two-layer model in 
the near wall region. 


Airfoil with Rime Ice 

Figures 16 and 17 show the results obtained by using 
the WIND code with the S-A turbulence model for the 
GLC305 airfoil with the 212 rime ice. These figures show 
CFD to predict accurately the lift, drag, and pressure 
coefficients until near stall. This shows CFD predictions 
of airfoils with rime ice, where separated regions are small, 
to be quite reliable. Similar to the situation with the clean 
airfoil, the stall angle is slightly under predicted. Again, 
the reason is that the flow becomes unsteady near stall, and 
only steady-state solutions were generated in this study. 

Prediction of the Velocity Field 

Figures 18 to 20 show the predicted contours of the x- 
component velocity magnitude for the GLC305 airfoil with 
944 glaze ice at AOA equal to 0°, 4°, and 6°. Also, shown 
are the experimentally measured values reported in Ref. 
12. Comparing the CFD results with the measured ones 
shows that CFD incorrectly predicts the size of the 
separated region downstream of the hom on airfoil’s 
suction side. This caused the surface pressure to be 
shifted, which in turn caused lift and drag to be under 
predicted. 

Figures 21 and 22 show the CFD and measured 
contours of the x-component velocity magnitude for the 
GLC305 airfoil with 212 rime ice at AOA equal to 6° and 
8°. These figures show the separated region to be 
predicted correctly though there are still discrepancies 
between the CFD and the measured results. Since the lift, 
drag, and pressure coefficients were predicted well by 
CFD, this indicates that predicting the separated region 
correctly is paramount. 

5. SUMMARY 

This study showed that if there are no large separated 
regions (e.g., the 212 ice), then CFD based on the one- 
equation Spalart-Allmaras turbulence model can predict 
accurately the lift, drag, and pressure coefficients. Thus, 
CFD can predict airfoils with rime ice quite adequately for 
angles of attack up to near stall. For airfoils with glaze ice 
(e.g., 944 ice), where the horns produce large separated 
regions about the airfoil’s leading edge, CFD predictions 
are much less satisfactory even at low angles of attack. 
For airfoils with glaze ice, this study showed that even the 
v 2 -f and the differential Reynolds stress models do not 
provide better results than the simple S-A model, which 
was found to provide the best results. However, more 
study is needed for the v 2 -f and the differential Reynolds 
stress models in which the near-wall treatment is not the 
two-layer model of Chen and Patel. Comparing the 
predicted x-component velocity magnitude with the 
measured ones show CFD to over predict the size of the 
separated region induced by the horns, and hence 
incorrectly predicting lift, drag, and pressure coefficients. 



This study also showed that for glaze ice, some 
smoothing of the ice shape is acceptable, which makes 
grid generation easier. 

Finally, it is noted that WIND and Fluent provided 
nearly identical results for lift, drag, and pressure 
coefficients when the following were the same: grid, 
turbulence model (S-A), and similar order of accuracy for 
the solution algorithms. 
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Fig. 1. NLF0414 airfoil with 623 ice. (a) Grid lines 
clustered next to solid boundaries propagate into domain 
along block boundaries, (b) Grid with wrap-around grid. 



Fig. 8. Computed and measured lift coefficient for clean 
GLC305 airfoil (M = 0.12, P = 20.5 psi, Re - 3.5 x 10 6 ). 
CFD: WIND, S-A model. 
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Fig. 9. Computed and measured pressure coefficient for clean 
GLC305 airfoil (M - 0.12, P = 20.5 psi, Re = 3.5 x 10 6 ). Top: 
AOA = 4°. Bottom: AOA = 6°. CFD: Wind, S-A model. 



Fig. 10. Computed and measured lift coefficient for 
GLC305 airfoil with 944 ice (M = 0.12, P = 20.5 psi, Re = 
3.5 x 10 6 ). CFD: WIND, S-A model. 
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Fig. 1 1 . Computed and measured pressure coefficient for GLC305 
airfoil with 944 ice (M = 0.12, P - 20.5 psi, Re = 3.5 x 10 6 ). Top: 
AOA = 4°. Bottom: AOA = 6°. CFD: Wind, S-A model. 
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Fig. 13. Effects of grid resolution and smoothing 944 ice 
shape and grid (M = 0.12, P = 20.5 psi, Re = 3.5 x 10 6 ). 
CFD: WIND, S-A model. 
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Fig. 15. Computed and measured lift and drag coefficient for 
GLC305 airfoil with 944 ice at AOA = 4° and 6° (M = 0.12, P 
= 37.0 psi. Re = 6.0 x 10 6 ). 








Fig. 16. Computed and measured lift and drag 
coefficients for GLC305 airfoil with 212 ice (M - 0.12, P 
= 20.5 psi, Re - 3.5 x 10 6 ). CFD: WIND, S-A model. 
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Fig. 17. Computed and measured pressure coefficient for 
GLC305 airfoil with 212 ice as a function normalized distance 
along the chord (M = 0.12, P = 20.5 psi, Re - 3.5 x 10 6 ). Top: 
AOA - 4°. Bottom: AOA = 6°. CFD: Wind, S-A model. 
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Fig. 18. Normalized x-component velocity magnitude for 
GLC305 airfoil with 944 ice at AOA = 0° (M = 0.12, P - 
20.5 psi, Re = 3.5 x 10 6 ). Top: CFD with WIND and S-A 
model. Bottom: measurement (Ref. 12). 


Fig. 19. Nonnalized x-component velocity magnitude for 
GLC305 airfoil with 944 ice at AOA = 4° ( M = 0.12, P = 
20.5 psi, Re = 3.5 x 10 6 ). Top: CFD with WIND and S-A 
model. Bottom: measurement (Ref. 12). 








Fig. 20. Normalized x-component velocity magnitude for 
GLC305 airfoil with 944 ice at AOA = 6° (M = 0.12, P = 
20.5 psi, Re = 3.5 x 10 6 ). Top: CFD with WIND and S-A 
model. Bottom: measurement (Ref. 12). 


Fig. 22. Normalized x-component velocity magnitude for 
GLC305 airfoil with 212 ice at AOA = 8° (M = 0.12, P = 20.5 
psi, Re - 3.5 x 10 6 ). Top: CFD with WIND and S-A model. 
Bottom: measurement (Ref. 12). 
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Fig. 21. Normalized x-component velocity magnitude for 
GLC305 airfoil with 212 ice at AOA = 6° (M = 0.12, P = 
20.5 psi, Re — 3.5 x 10 6 ). Top: CFD with WIND and S-A 
model. Bottom: measurement (Ref. 12). 




